Infinite-randomness critical point in the two-dimensional disordered contact process 



Thomas Vojta, Adam Farquhar, and Jason Mast 
Department of Physics, Missouri University of Science and Technology, Rolla, MO 65409, USA 

(Dated: January 13, 2009) 

We study the nonequihbrium phase transition in the two-dimensional contact process on a ran- 
domly diluted lattice by means of large-scale Monte-Carlo simulations for times up to 10^° and 
system sizes up to 8000 x 8000 sites. Our data provide strong evidence for the transition being 
controlled by an exotic infinite-randomness critical point with activated (exponential) dynamical 
scaling. We calculate the critical exponents of the transition and find them to be universal, i.e., 
independent of disorder strength. The Griffiths region between the clean and the dirty critical 
points exhibits power-law dynamical scaling with continuously varying exponents. We discuss the 
generality of our findings and relate them to a broader theory of rare region effects at phase tran- 
sitions with quenched disorder. Our results are of importance beyond absorbing state transitions 
because according to a strong-disorder renormalization group analysis, our transition belongs to the 
universality class of the two-dimensional random transverse-field Ising model. 
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I. INTRODUCTION 

Many-particle systems far from equilibrium can un- 
dergo phase transitions if external parameters are varied. 
These transitions, which separate different nonequihb- 
rium steady states, are characterized by large-scale fluc- 
tuations and collective behavior over large distances and 
long times very similar to the behavior at equilibrium 
critical points. Examples of nonequilibrium phase tran- 
sitions occur in a wide variety of systems ranging from 
surface chemical reactions and growing interfaces to traf- 
fic jams and to the spreading of epidemics in biology. 
Reviews of some of these transitions can be found, e.g., 
in Refs. [1, 2, 3, 4, 5, 6, 7] 

In recent years, a large effort has been directed towards 
classifying possible nonequilibrium phase transitions [5]. 
A particularly well-studied type of transitions separates 
active, fluctuating steady states from inactive absorbing 
states where fluctuations cease entirely. The generic uni- 
versality class for these absorbing state transitions is di- 
rected percolation (DP) [8]. According to a conjecture 
by Janssen and Grassberger [9, 10], all absorbing state 
transitions with a scalar order parameter, short-range in- 
teractions, and no extra symmetries or conservation laws 
belong to this class. Examples include the transitions in 
the contact process [11], catalytic reactions [12], interface 
growth [13], or Pomeau's conjecture regarding turbulence 
[14]. In the presence of conservation laws or additional 
symmetries, other universality classes such as the parity 
conserving class or the Z2-symmetric directed percola- 
tion (DP2) class can occur (see Ref. [4] and references 
therein) . 

The DP universality class is ubiquitous in theory and 
computer simulations, but clearcut experimental verifica- 
tions were lacking for a long time [15]. Partial evidence 
was manifest in the spatio-temporal intermittency in fcr- 
rofluidic spikes [16]. Only very recently, a full verification 
(the only one, to the best of our knowledge) was found 
in the transition between two turbulent states in a liquid 



crystal [17]. One possible cause for the surprising rarity 
of DP scaling in experiments are impurities, defects, or 
other kinds of disorder that are present in all realistic 
systems. 

For this reason, the effects of quenched spatial disor- 
der on the DP universality class have been studied for 
many years. However, a consistent picture has been very 
slow to emerge. According to the general Harris criterion 
[18] (applied to the DP imiversality class by Kinzel [19] 
and Noest [20]), a clean critical point is (perturbatively) 
stable against weak spatial disorder, if the spatial corre- 
lation length critical exponent fulfills the inequality 
dv± > 2; where d is the spatial dimensionality. The cor- 
relation length exponents of the DP universality class are 
jy^ ~ 1.097 in one space dimension, 0.73 in two dimen- 
sions, and 0.58 in three dimensions [4], therefore, the Har- 
ris criterion is violated in all dimensions d < A. The insta- 
bility of the DP critical behavior against spatial disorder 
was confirmed by a field-theoretic renormalization group 
calculation [21]. This study did not produce a new criti- 
cal fixed point but runaway flow towards large disorder, 
suggesting unconventional behavior [55]. Early Monte- 
Carlo simulations in two space dimensions [23, 24] found 
logarithmically slow dynamics in violation of power-law 
scaling. Moreover, in analogy with Griffiths singularities 
[25], very slow dynamics was found in an entire parame- 
ter region near the transition [26, 27, 28, 29]. 

A crucial step towards resolving this puzzling sit- 
uation was taken by Hooyberghs et al. [30, 31] who 
used the Hamiltonian formalism [32] to map the one- 
dimensional disordered contact process onto a random 
quantum spin chain. They then applied a Ma-Dasgupta- 
Hu strong-disorder renormalization group [33, 34] and 
showed that the transition is controlled by an exotic 
infinite-randomness fixed point in the universality class 
of the random transverse-field Ising model [35, 36], at 
least for sufficiently strong disorder. This type of criti- 
cal point displays ultraslow activated rather than power- 
law dynamical scaling. For weaker disorder, Hooyberghs 
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et al. relied on numerical simulations and predicted 
non-universal continuously varying exponents, with ei- 
ther power-law or exponential dynamical scaling. Using 
large-scale Monte-Carlo simulations Vojta and Dickison 
[37] recently confirmed the infinite-randomness scenario. 
Moreover, their critical exponent values agreed with the 
predictions of the strong-disorder rcnormalization group 
(which can be solved exactly in one dimension) , and they 
were universal, i.e., independent of disorder strength. 

In higher dimensions, a strong-disorder renormaliza- 
tion group can still be applied, but, in contrast to one di- 
mension, it cannot be solved analytically. By implement- 
ing the rcnormalization group numerically, Motrunich 
ct al. [38] demonstrated the existence of an infinite- 
randomness fixed point in two dimensions (and at least 
the flow towards larger disorder in three dimensions). 
However, reliable estimates for the critical exponent val- 
ues have been hard to come by due to the small sizes 
available in the quantum simulations. 

In this paper, we report the results of large-scale 
Monte-Carlo simulations of the contact process on two- 
dimensional randomly diluted lattices. The purpose of 
this work is twofold. First, we intend to determine be- 
yond doubt the character and universality of the critical 
behavior; is it of conventional power-law type, infinite- 
randomness type or an even more exotic kind? Sec- 
ond, we intend to perform a careful data analysis that 
allows us to compute reasonably accurate exponent val- 
ues. These will be important beyond the disordered con- 
tact process and apply to all systems controlled by the 
same infinite-randomness fixed point including the two- 
dimensional random transverse-field Ising magnet [38] . 

The paper is organized as follows. In section II, we 
introduce our model, the contact process on a randomly 
diluted lattice. We then contrast the scaling theories for 
conventional critical points and infinite-randomness crit- 
ical points. We also summarize the predictions for the 
Griffiths region. In section III, we present our simula- 
tion method and various numerical results together with 
a comparison to theory. We conclude in section IV by 
summarizing our findings, discussing three space dimen- 
sions, and relating our results to a broader classification 
of phase transitions with quenched disorder [39]. 



II. MODEL AND PHASE TRANSITION 
SCENARIOS 

A. Contact process on a diluted lattice 

The contact process [11], a prototypical system in the 
DP universality class, can be interpreted as a simple 
model for the spreading of a disease. The clean contact 
process is defined on a d-dimensional hypercubic lattice. 
Each lattice site r can be active (infected) or inactive 
(healthy). In the course of the time evolution, active sites 
can infect their neighbors or they can spontaneously heal 
(become inactive). Specifically, the dynamics of the con- 



tact process is given by a continuous-time Markov pro- 
cess during which active sites become inactive at a rate fi 
while inactive site become active at a rate An/ (2d). Here, 
n is the number of active nearest neighbor sites. The in- 
fection rate A and the healing rate /i (which can be set to 
unity without loss of generality) are external parameters. 
Their ratio controls the behavior of the system. 

For A /i, healing dominates, and the absorbing state 
without any active sites is the only steady state (inactive 
phase). For sufficiently large infection rate A, there is a 
steady state with a nonzero density of active sites (active 
phase) . These two phases are separated by a nonequilib- 
rium phase transition in the DP universality class at a 
critical infection rate A|^. 

We now introduce quenched spatial disorder into the 
contact process by randomly diluting the underlying lat- 
tice. Specifically, we randomly remove each lattice site 
with probability p [56] . As long as the impurity concen- 
tration p remains below the percolation threshold Pc, the 
lattice has an infinite connected cluster of sites which can 
support the active phase. Conversely, at dilutions above 
Pc, the lattice is decomposed into disconnected finite-size 
clusters. Because the infection eventually dies out on any 
finite cluster, there is no active phase for p > pc- Conse- 
quently, the contact process on a randomly diluted lattice 
has two different nonequilibrium transitions (see also Fig. 
1), the generic transition for p < pc which is driven by 
the dynamic fluctuations of the contact process and a 
percolation transition ai p = Pc which is driven by the 
lattice geometry [41]. The two transitions are separated 
by a multicritical point which was studied in Rcf. [42] . 

The basic observable in the contact process is the av- 
erage density of active sites at time <, 

/'W = ;^E("rW) (1) 

r 

where nr(t) = 1 if the site r is active at time t and 
rir(t) = if it is inactive. L is the linear system size, 
and (...) denotes the average over all realizations of the 
Markov process. The longtime limit of this density (i.e., 
the steady state density) 

Pstat = lim p{t) (2) 

t — 'OO 

is the order parameter of the nonequilibrium phase tran- 
sition. 



B. Conventional power-law scaling 

In this subsection, we briefly summarize the scaling 
theory for absorbing state transitions controlled by con- 
ventional fixed points with power-law dynamical scaling 
(see, e.g., Ref. [4]), using the clean contact process as an 
example. 

When the critical point Ac is approached from the ac- 
tive phase, the order parameter (steady state density) 
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vanishes following the power law 

fttat (A - ~ (3) 

where A = (A — Ac)/ Ac is the dimensionless distance from 
the critical point, and /3 is the order parameter critical 
exponent. In addition to the average density, we also 
need to characterize the length and time scales of the 
density fluctuations. When approaching the transition, 
the (spatial) correlation length diverges as 

u - \^r^ • (4) 

The correlation time ^|| behaves like a power of the cor- 
relation length, 

eii-ei, (5) 

i.e., the dynamical scaling is of power-law form. The 
three critical exponents /3, I'x and z completely char- 
acterize the directed percolation universality class. The 
above relations allow us to write down the finite-size scal- 
ing form of the density as a function of A, t, and L, 

p{^,t,L) ^b^^'"^ pi^b-^'"^ ,tb\Lb) , (6) 

where b is an arbitrary dimensionless scale factor. 

Two interesting observables can be studied if the time 
evolution starts from a single active site in an other- 
wise inactive lattice. The survival probability Ps{t) is 
the probability that an active cluster survives at time 
t when starting from such a single-site seed at time 0. 
For directed percolation, the survival probability scales 
exactly like the density [57], 

P,{A,t,L) ^ bl^^^^P.iAb-^/"^ ,tt\Lb) . (7) 

The pair connectedness function C(r,t) = {n.r{t) no{0)) 
describes the probability that site r is active at time t 
when starting from an initial condition with a single ac- 
tive site at r = and time 0. Because C involves a 
product of two densities, its scale dimension is 2(3 /i'±, 
and the full finite-size scaling form reads [58] 

C{A,r,t,L) = b^^^''^C{Ab-^/''^,rb,tb'-,Lb) . (8) 

The total number of particles N when starting from a 
single seed site can be obtained by integrating the pair 
connectedness C over all space. This leads to the scaling 
form 

N{A,t,L) = b'^'^'''^-'^N{Ab-~^/''^,tb\Lb) . (9) 

Finally, the mean-square radius R of the active cluster, 
being a length, scales as 

i?(A, i, L) = b-^RiAb-^/"^ , ^6^ Lb) . (10) 

At the critical point, A = 0, and in the thermody- 
namic limit, L — > cxj, the above scaling relations lead 
to the following predictions for the time dependencies 



of observables: The density and the survival probability 
asymptotically decay like 

p{t)^r\ Ps{t)^r' (11) 

with 5 = P/{i^±z). In contrast, the radius and the num- 
ber of particles in a cluster starting from a single seed 
site increase like 

i^(^)-^l/^ 7V(t)-i® (12) 

where Q = d/z — 2(3 /{v^z) is the so-called critical initial 
slip exponent. 

C. Activated scaling at an infinite- randomness 
fixed point 

In this subsection we summarize the scaling theory for 
an infinite-randomness fixed point with activated scal- 
ing, as has been found in the one-dimensional disordered 
contact process [30, 37, 39]. 

At an infinite-randomness fixed point, the dynamics 
is ultraslow. The power-law dynamical scaling (5) gets 
replaced by activated scaling 

MeiiAo)~Ct, (13) 

characterized by a new exponent i}} called the tunnel- 
ing exponent. (This name stems from the random 
transverse-field Ising model where this type of scaling 
was first found.) Here to is a nonuniversal microscopic 
time scale. The exponential relation between time and 
length scales implies that the dynamical exponent z is 
formally infinite. In contrast, the static scaling behavior 
remains of power law type. 

Another important characteristic of an infinite- 
randomness fixed point is that the probability distribu- 
tions of observables become extremely broad. Therefore, 
the average and typical value of an observable do not 
necessarily agree because averages may be dominated by 
rare events. Nonetheless, the scaling form of the aver- 
age density at the infinite-randomness critical point is 
obtained by simply replacing the power-law scaling com- 
bination tV' by the activated combination h\{t/to)b^ in 
the argument of the scaling function: 

p(A,ln(t/to),L) =5'3/-XA6-i/'^\ln(t/to)5'^,L6) . 

(14) 

Analogously, the scaling forms of the average survival 
probability, the average number of sites in a cluster, and 
the mean-square cluster radius starting from a single site 
are 

Ps{AM{t/t^),L) = b''/''^P,{Ab-'/-'\Ht/to)b^,Lb) 

(15) 

iV(A, ln(t/to), L) = b^'^/^^-'^NiAb-^^"^ . \n{t/to)b'^ , Lb). 

(16) 
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i?(A,ln(t/to),i) = b-^ RiAb-^/"^ ,lii{t/to)b'f' , Lb). 

(17) 

These activated scaling forms lead to logarithmic time 
dependencies at the critical point (in the thermodynamic 
limit). The average density and the survival probability 
asymptotically decay like 

p{t) ^ Mt/to)]-', Ps{t) ^ Mt/to)]-' (18) 

with S = (3 /(y±.4') while the radius and the average num- 
ber of particles in a cluster starting from a single seed 
site increase like 

R{t) ~ Mt/t^)Y'^, Nit) ^ Mt/to)f (19) 

with e ^d/ij- 2l3/{v±ij). 

In contrast to one dimension, where the critical ex- 
ponents can be calculated exactly within the strong- 
disorder renormalization group, they need to be found 
numerically in higher dimensions. 



D. Griffiths singularities 

Quenched spatial disorder does not only destabilize the 
DP critical point, it also leads to interesting singularities, 
the so-called Griffiths singularities [25], in an entire pa- 
rameter region around the transition. For a recent review 
of this topic at thermal, quantum, and nonequilibrium 
transitions, see Ref. [39]. In the context of the contact 
process on a diluted lattice, Griffiths singularities can be 
understood as follows. 

Because lattice dilution reduces the tendency towards 
the active phase, the critical infection rate of the diluted 
system is higher than that of the clean system, Ac > A°. 
For infection rates A with Ac > A > A°, the system is 
globally in the inactive phase, i.e., it eventually decays 
into the absorbing state. However, in the thermodynamic 
limit, one can find arbitrarily large spatial regions devoid 
of impurities. These rare regions are locally in the active 
phase. They cannot support a non-zero steady state den- 
sity because they are of finite size, but their decay is very 
slow because it requires a rare, exceptionally large den- 
sity fluctuation. 

As a result, the inactive phase of the contact process 
on a diluted lattice can be divided into two regions. For 
infection rates below the clean critical point, A < Aj?, 
the behavior is conventional. The system approaches the 
absorbing state exponentially fast in time. The decay 
time increases with A and diverges as |A — Acl"^"^-^ where 
z and i'± are the exponents of the clean critical point 
[37, 45]. Inside the Grifhths region, Ac > A > A^l, one 
has to estimate the rare region contribution to the time 
evolution of the density [20, 26]. The probability w for 
finding a rare region of linear size Lr devoid of impurities 
is (up to pre-exponcntial factors) given by 

w{Lr) - cxp(-pL^) (20) 



with p = — ln(l — p). To exponential accuracy, the rare 
region contribution to the density can be written as 

p(t) - J dLr Lf w{Lr)exp[~t/T{Lr)] (21) 

where T{Lr) is the decay time of a rare region of size 
Lr- Let us first discuss the behavior at the clean critical 
point, Ac, i.e., at the boundary between the conventional 
inactive phase and the Griffiths region. Here, the decay 
time of a single, impurity-free rare region of size Lr scales 
as t{Lj.) ^ L^ with z the clean critical exponent [see (6)]. 
Using the saddle point method to evaluate the integral 
(21), we find the leading long-time decay of the density 
to be given by a stretched exponential, 

\np{t) - t'^/(d+^) , (22) 

rather than a simple exponential decay as for A < Aj!. 

Within the Griffiths region, A^ < A < Ac, the decay 
time of a single rare region depends exponentially on its 
volume, 

T{Lr) - exp(aLf) (23) 

because a coordinated fluctuation of the entire rare region 
is required to take it to the absorbing state [20, 26, 46]. 
The nonuniversal prefactor a vanishes at the clean critical 
point Ac and increases with A. Close to A^, it behaves 
as a ^J'' (A — Ac)'*''^ with vj_ the clean critical 
exponent. Repeating the saddle point analysis of the 
integral (21) for this case, we obtain a power-law decay 
of the density 

pit) ~ t-P'" ^ t-"^/'' (24) 

where z' — da/p is a customarily used nonuniversal dy- 
namical exponent in the Griffiths region. Its behav- 
ior close to the dirty critical point Ac can be obtained 
within the strong disorder renormalization group method 
[31, 36, 38]. When approaching the phase transition, z' 
diverges as 

z' ~ |A- Ac]-^''^ (25) 

where ip and vj^ arc the exponents of the dirty critical 
point. 

Let us emphasize that strong power-law Griffiths sin- 
gularities such as (24) usually occur in connection with 
infinite-randomness critical points, while conventional 
critical points display exponentially weak Griffiths effects 
[39]. Thus, the character of the Griffiths singularities can 
be used to identify the correct scaling scenario. 

III. MONTE-CARLO SIMULATIONS 

A. Metliod and overview 

We now turn to the main objective of our work, large- 
scale Monte-Carlo simulations of the contact process on 
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a randomly diluted square lattice. There are several effi- 
cient computational implementations of the contact pro- 
cess, all equivalent with respect to the universal behavior 
at the phase transition. We followed the algorithm de- 
scribed, e.g., by Dickman [47]. The simulation starts at 
time t = from some configuration of active and inactive 
sites. Each event consists of randomly selecting an active 
site r from a list of all Na active sites, selecting a process: 
creation with probability A/(l -I- A) or annihilation with 
probability 1/(1 -I- A) and, for creation, selecting one of 
the four neighboring sites of r. The creation succeeds, if 
this neighbor is empty (and not an impurity site). The 
time increment associated with this event is l/Na- 

Using this method, we investigated systems with sizes 
of up to 8000 X 8000 sites and impurity (vacancy) concen- 
trations p = Q, 0.1, 0.2, 0.3, and p = Pc ^ 0.407254. To 
explore the ultraslow dynamics predicted in the infinite- 
randomness scaling scenario, we simulated very long 
times up to i = 10^" which is, to the best of our knowl- 
edge, significantly longer than all previous simulations 
of the disordered two-dimensional contact process. In 
all cases, we averaged over many disorder realizations, 
details will be mentioned below for each specific simula- 
tion. The total numerical effort was approximately 40000 
CPU-days on the Pegasus Cluster at Missouri S&T. 

To identify the critical point for each parameter set, 
and to determine the critical behavior, we carried out two 
types of simulations, (i) runs starting from a completely 
active lattice during which we monitor the time evolution 
of the density p{t), and (ii) runs starting from a single 
active site in an otherwise inactive lattice during which 
we monitor the survival probability Ps(t) and the size 
Ns{t) of the active cluster. Figure 1 gives an overview 
of the phase diagram resulting from these simulations. 
As expected, the critical infection rate Ac increases with 
increasing impurity concentration. In the following sub- 
sections we discuss the behavior in the vicinity of the 
phase transition in more detail. 
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FIG. 1: (Color online) Phase diagram of the contact process 
on a site-diluted square lattice (inverse critical infection rate 
A^^ vs impurity concentration p). MCP marks the multicrit- 
ical point. The black dots show the actual simulation results, 
the lines are guides to the eye. 
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FIG. 2: (Color online) Time evolution of the average density 
p{t) for the contact process on an undiluted (p = 0) square 
lattice, starting from a completely active lattice. 



B. Clean 2d contact process 

We first performed a number of simulations for clean 
undiluted lattices {p = 0), mainly to test our numerical 
implementation of the contact process. Fig. 2 shows the 
density p{t) for runs starting from a completely active 
lattice. The data are averages over 60 runs of a system 
of linear size L = 4000 except for the two infection rates 
closest to the transition, A = 1.64874 and 1.64875 for 
which we performed 42 runs of a system of size L = 8000. 
From these simulations, we estimate the critical infection 
rate to be Aj? = 1.64874(2) where the number in brackets 
gives the error of the last digit. The critical exponent S 
can be determined from a power-law fit of p{t); we find 
S = 0.4526(7). 

We then perform a scaling analysis of p{t) [based on 
(6)] on the inactive side of the transition using 12 dif- 
ferent infection rates between 1.60 and A^ = 1.64874. 



This allows us to find the exponent combination zu± = 
1.290(4). Finally, we carried out several runs right at the 
critical point but with system sizes varying from L = 25 
to L = 8000. A finite-size scaling analysis [again based 
on (6)] gives the dynamical exponent z = 1.757(8). All 
other critical exponents can now be calculated from the 
scaling relations (6), (7), and (9). We find = 0.734(6), 
/? = 0.584(3), and 9 = 0.233(6). 

To verify the results we also carried out simulations 
starting from a single active seed site. Using 500,000 
runs with a maximum time of 10^, we arrived at the 
same value A° = 1.64874(2) for the critical infection rate. 
We measured the survival probability Ps (t) , the number 
of active sites Ns{t) and the mean square radius of the 
active cluster R{t). The resulting estimates of the expo- 
nents S, 0, and z confirmed the values quoted above. 

Our results are in excellent agreement with, and of 
comparable accuracy to, other large-scale simulation of 
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FIG. 3: (Color online) Time evolution of the survival prob- 
ability Ps(t) for impurity concentration p = 0.2, starting 
from a single seed site. The clean critical infection rate 
Ac = 1.64874 and the actual critical point Ac = 2.10750 are 
specially marked. 
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the DP universality class in two dimensions (see, e.g., 
Ref. [47] and references therein). 
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FIG. 4: (Color online) InPs(t) and \nNa{t) versus In Int for 
several infection rates close to the critical point (at least 5000 
disorder realizations with 128 trials each, p = 0.2). 



C. Diluted 2d contact process: Overview 

We start our discussion of the diluted contact process 
by showing in Fig. 3 an overview of the survival prob- 
ability Ps{t) for a system with impurity concentration 
p = 0.2, covering the A range from the conventional in- 
active phase, A < A°, all the way to the active phase, 
A > Ac. The data represent averages over 5000 disorder 
configurations, with 128 runs for each disorder configu- 
ration starting from random seed sites. The system size, 
L = 2000, is chosen such that the active cluster stays 
much smaller then the sample over the time interval of 
the simulation (thus eliminating finite-size effects). 

For infection rates below and at the clean critical point 
A° = 1.64874, the density decay is very fast, clearly faster 
than a power law. Above A°, the decay becomes slower 
and asymptotically appears to follow a power law with 
an exponent that varies continuously with A. For even 
larger infection rates, the decay seems to be slower than a 
power law. At the largest infection rates, the system will 
ultimately reach a nonzero steady-state survival proba- 
bility, i.e., the it is in the active phase. 

Let us emphasize that the behavior in Fig. 3, in par- 
ticular, the nonuniversal power-law decay over a range of 
infection rates, already provides support for the infinite- 
randomness scenario of subsection II C while it disagrees 
with the conventional scenario (subsection II B) of power- 
law behavior restricted to the critical point and exponen- 
tially weak Griffiths effects. 



D. Finding the critical point 

A very efficient method for identifying the critical point 
in the clean contact process and other clean reaction- 
diffusion systems is to look for the critical power-law time 
dependencies of various observables. To use the same 
idea in our diluted contact process simulations, we plot 
InPs(t), \nNs{t) and In /?(<:) versus In Int. In such plots, 
the anticipated logarithmic laws (18) and (19) correspond 
to straight lines. Figure 4 shows these plots for several 
infection rates close to the critical point. The parameters 
are as in subsection IIIC. At the first glance, the data 
at the appropriate critical infection rate appear to follow 
the expected logarithmic behavior over several orders of 
magnitude in time t. However, a closer inspection reveals 
a serious problem. The survival probability Pg suggests a 
critical infection rate Ac in the range 2.10500 to 2.10625 
while the cluster size Ng appears subcritical even at the 
larger value 2.10750. Such a discrepancy suggests strong 
corrections to scaling. 

The origin of these corrections can be easily under- 
stood by rewriting Pg{t) - [ln(i/to)]"^ = [ln(t) - 
ln(io)]~''. The microscopic time scale thus provides a 
correction to scaling; and since we cover only a mod- 
erate range in ln(t), it strongly influences the results. 
We emphasize that this problem is intrinsic to the case 
of activated scaling (and thus unavoidable) because log- 
arithms, in contrast to power laws, are not scale free. 
Identifying straight lines in plots like Fig. 4 is thus not 
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FIG. 5; (Color online) lnNs{t) versus \nPs{t) and for several 
infection rates close to the critical point. The straight line 
is a power-law fit of the asymptotic part of the A = 2.10750 
curve, (p = 0.2, 5000 - 20,000 disorder reahzations with 128 
trials each. 

a reliable tool for finding the critical point, even more so 
because simulations in ID (where the analysis is much 
easier because the critical exponents are known exactly) 
have given sizable to ~ 100 . . . 500 for typical parameter 
values [37]. 

To circumvent this problem, we devised a method for 
finding the critical point without needing a value for the 
microscopic scale to . It is based on the observation that in 
the infinite-randomness scenario, to has the same value in 
the scaling forms of all observables (because it is related 
to the basic energy scale of the strong-disorder renor- 
malization group). Thus, if we plot Ns{t) versus Ps{t), 
the critical point corresponds to power-law behavior pro- 
vided all other corrections to scaling are weak. The same 
is true for other combinations of observables. 

Figure 5 shows a plot of Ns (t) versus Pg (t) for several 
infection rates close to the phase transition. The figure 
shows that the data for A = 2.1075 follow the expected 
power law for about a decade in P, (10~^ < Pg < 10"^, 
corresponding to 4 decades of time, from about 10'^ to 
10*), while the data at higher and lower A curve away 
as expected. We therefore identify Ac = 2.1075(1) as the 
critical infection rate for dilution p = 0.2. This value 
agrees well with an estimate based on Dickman's [23] 
heuristic criteria of Ac being the smallest A supporting 
asymptotic growth of Ns{t). 

Figure 5 also reveals two remarkable features of the 
problem, (i) The asymptotic critical region is extremely 
narrow, |A| = |A — Ad /Ac ~ 10""^, at least for these 
two observables. This follows from fact that the data 
for A outside the small interval [2.1070, 2.1080] curve 
away from the critical curve well before it reaches the 
asymptotic regime, (ii) The crossover to the asymptotic 
behavior occurs very late since P^ ~ 0.01 corresponds 
to a crossover time oi t^ ~ 10^, implying that very long 
simulations are necessary to extract the true critical be- 
havior. Moreover, the mean-square radius of a surviving 



active cluster at the crossover time is R{tx) ~ 55, cor- 
responding to an overall diameter of about 200. This 
means simulations with linear system sizes below about 
200 will never reach the asymptotic behavior. We con- 
firmed this observation by carrying out simulations with 
smaller sizes. 

Having determined the critical point, we can now ob- 
tain a rough estimate of the microscopic time scale to by 
replotting lnPs{t) and lnNs{t) as functions of\n\n(t/to) 
with varying to- At the correct value for to-, the critical 
curves (A = 2.1075) of all observables must be straight 
lines. By performing this analysis for Ps{t), Ns{t), and 
p{t) we find Into = 5-5 ± 1.0 for dilution p = 0.2. The 
results for all the observables agree within the error bars, 
providing an additional consistency check for our method 
and the underlying infinite-randomness scenario. 

E. Critical behavior 

According to subsection II C, the critical behavior at 
the infinite-randomness fixed point can be characterized 
by three independent exponents, e.g., v^, and ^j. The 
two "static" exponents (3 and i^^ can in principle be de- 
termined without reference to to, the value of the tunnel- 
ing exponent ip depends on the above estimate for to- 

The combination (3/i'± (the scale dimension of the or- 
der parameter) can be obtained directly from the critical 
curve in Fig. 5. By combining cqs. (18) and (19), we 
find Ns ~ Ps~^^^^ ~ p2ii-^±/i3)_ ^ ^j^^ asymptotic 

part of the critical curve gives (3/i^_l = 0.96(2). Here, 
the error is almost entirely due to the uncertainty in lo- 
cating the critical infection rate Ac, the statistical error 
is much smaller. We obtained independent estimates for 
P/v^ from analyzing plots of Ng vs. R and Pg vs. R in 
a similar fashion. The results agree with the above value 
within their error bars. 

To obtain the tunneling exponent ■0, we now fit the 
critical survival probability Pg{t) to [ln(t/to)]~^ (for 
times up to t = lO'^) using our estimate ln<o = 5.5 ± 1.0. 
We find 5 = 1.9(2) with the error mainly coming from the 
uncertainties in Ac and to- An analogous analysis of the 
density of active sites p{t) in a simulation starting from 
a fully infected lattice (L = 8000, times up to t = 10^) 
gives exactly the same value (see Fig. 6). From the time 
evolution of Ns{t), we also find the estimate 8 = 0.15(3) 
for the initial slip exponent. Using tjj = P/{h'±6) gives 
the tunneling exponent, tp = 0.51(6). Estimates obtained 
from Ng{t) and directly from fitting R{t) to [ln(t/to)]^/''' 
agree within the error bars with this value. It must be 
emphasized that the value of tp sensitively depends on 
the correct identification of to- A naive fit of Pg{t) to a 
power of ln(t) would give a value of about 3.0 for S and 
thus a much smaller value of about 0.33 for ip. 

Finding the final missing exponent in our triple (/?, i'^ , 
ijj) requires the analysis of off-critical simulation runs. 
Here, the extreme narrowness of the asymptotic criti- 
cal region causes serious problems, because wc have only 
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FIG. 6: (Color online) In p(t) versus lnln(f/to) and for several 
infection rates at and below to the critical point. The dashed 
line represents 2/3 of the critical p. Inset: Time at which 
the density falls below 2/3 of the critical value as a function 
of the distance from criticality. The solid line is a power-law 
fit of ln(fa:/to) vs A. (p = 0.2, 13 disorder realizations with 
L = 8000 for A > 2.1, 250 realizations with L = 2000 for 
A < 2.1) 



FIG. 7: (Color online) Inp(f) vs. InA(f) at the crossover 
time t — tx for the density itself and \nPs{t) vs. InA(t) at 
the crossover time t = tx for the particle number Ns- The 
dashed lines are power-law fits to the expected behavior p ~ 
~ lAI''. 



F. Universality 



three A values from 2.1070 to Ac = 2.1075. Moreover, a 
glance at Fig. 5 shows that adding simulations runs at 
additional A closer to the critical point would only make 
sense if we also significantly increased the simulation time 
and/or the number of disorder realizations. Because this 
appears to be beyond our present computational capa- 
bilities, we determine an estimate for P by extrapolating 
the effective exponent obtained outside the asymptotic 
critical region. To this end, we analyze the time evo- 
lution of the density p(t) starting from a fully infected 
lattice show in Fig. 6 for various infection rates A below 
the critical point. For each A, we define a crossover time 
tx at which the density p{t) is exactly 2/3 of the critical 
density at the same time. 

This crossover time can be analyzed in two ways. Ac- 
cording to (14), it should behave as \n{tx/to) ^ lA]'"-^^ . 
The inset of Fig. 6 shows a fit to this form, giving 
iy±ip = 0.64(10). However, this value depends on our 
estimate for the microscopic scale to . To avoid relying on 
to, we plot |A(i:i,)| vs. the critical p{t) (at t = t^) m Fig. 
7. According to (14), the expected behavior is a power 
law p ~ lAj'^. The same analysis can also be performed 
by analyzing the crossover of Ns(t) and plotting the re- 
sulting A(ta;) vs. Ps{t)- Both data sets show significant 
curvature (i.e. deviations from the expected power law 
behavior) which is not surprising because the underlying 
A are not in the asymptotic region. By fitting the small- 
A part of both curves, we estimate /? = 1.15(15). A fit 
of A{tx) vs. Nsit) gives i/^ - /3 = 0.045(15). Thus, 
our final estimate for the correlation length exponent 
i>±^ = 1.20(15) in agreement with the Harris criterion 
> 2/d= 1. 



The universality of the infinite-randomness scenario for 
the critical point in the disordered contact process has 
been controversially discussed in the past [31, 37, 48]. 
The underlying strong-disorder renormalization group 
becomes exact only in the limit of broad disorder dis- 
tributions and can thus not directly address the fate of a 
weakly disordered system. However, Hoyos [49] recently 
showed that within this renormalization group scheme, 
the disorder always increases under renormalization, even 
if it is weak initially. Moreover, the perturbative renor- 
malization group of Janssen [21], which is controlled for 
weak disorder, displays runaway flow towards large disor- 
der strength supporting a universal scenario in which the 
infinite-randomness fixed point controls the transition for 
all bare disorder strength. 

To address this question in our simulations, we repeat 
the above analysis for impurity concentrations p = 0.1 
and 0.3 albeit with somewhat smaller numbers of disorder 
realizations. We determine the critical infection rates 
from plots of Ns{t) vs. Ps{t) analogous to Fig. 5 and find 
Ac = 1.8462(3) for p = 0.1 and Ac = 2.4722(2) for p = 
0.3. Figure 8 shows the critical Ng vs. Pg curves for p = 
0.1, 0.2, and 0.3. In the long-time (low- Ps) limit, all three 
curves follow power-laws, and the resulting values for the 
exponent P/i^±_ agree within their error bars, P/i'i^ = 
0.97(3) for p = 0.1, 13/va_ = 0.96(2) for p = 0.2, and 
l3/u± = 0.96(3) for p = 0.3. 

Having fixed the critical point for all p, we determine 
the microscopic time scale to as in subsection HID. Its 
value changes significantly with p, it is Into ~ 6.5 (for 
p = 0.1), 5.5 (for p = 0.2), and 4.0 (for p = 0.3). Power- 
law fits of Ps{t) and Ns{t) to ln(t/to) give estimates for 
the exponents 6 and Q. Within their error bars, the 
values for p = 0.1 and 0.3 agree with the respective values 
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FIG. 8: (Color online) \nNs{t) versus \nPs{t) at criticality 
for impurity concentrations p = 0.1, 0.2, and 0.3. The dashed 
lines are power-law fits of the asymptotic parts (18000 disor- 
der configurations for p = 0.2, 4000 each for p = 0.1 and 0.3, 
128 spreading trials for each configuration). 



for p = 0.2 discussed in subsection HIE. 

Our simulations thus provide numerical evidence for 
the critical behavior of the disordered 2D contact process 
being universal, i.e., being controlled by the same infinite- 
randomness fixed point for all disorder strength. Note 
that the bare disorder can be considered weak for p ~ 0.1 
as the difference of the critical infection rate from its 
clean value is small, (Ac — A°)/A° w 0.12. 



G. Griffiths region 

In addition to the critical point, we have also studied 
the Griffiths region A° < A < Ac in some detail. Fig- 
ure 9 shows a double-logarithmic plot of the density p(t) 
(starting from a fully infected lattice) for several infec- 
tion rates in this region. For all these A, the long-time 
behavior of the density follows the expected power law 
over several order of magnitude in time and/or density. 
The values of the (non-universal) dynamical exponent z' 
can be determined by fitting the long-time behavior to 
eq. (24). The inset of Fig. 9 shows that z' diverges with 
A approaching the critical point Ac = 2.1075. Fitting 
this divergence to the power law (25) expected within 
the infinite-randomness scenario, we find ipv^ ~ 0.61 in 
good agreement with the value extracted from the scaling 
analysis of the density in subsection III E. Performing the 
same analysis with the survival probability data shown 
in Fig. 3 gives analogous results. 



H. Spatial correlations 

To study the spatial correlations in the diluted con- 
tact process close to criticality, we calculate the average 
equal-time correlation function G = '^j.{nr{t)no{t)) 



FIG. 9; (Color online) Inp(t) versus Int for several infection 
rates in the Griffiths region, A2 < A < Ac (p = 0.2, 13 disorder 
realizations with L — 8000 for A > 2.1, 250 realizations with 
L = 2000 for A < 2.1 ). The dashed lines are power-law fits 
of the long-time behavior to (24). Inset: Resulting dynamical 
exponent z' as a function of A. 
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FIG. 10; (Color online) InG(r) vs. Inr for several infection 
rates A (p = 0.2, 13 disorder realizations with L = 8000). 
G(r) is calculated from the configurations at late times of the 
time evolution (t > 4 x 10*). 



from the late-time configurations arising in the simula- 
tions starting from a fully infected lattice. Scaling ar- 
guments analogous to those leading to (8) suggest the 
scaling form 

G(A,r,ln(t/to)) = b^^/''^G{Ab-^^''^ ,rbM{t/to)b'l') . 

(26) 

To find the true stationary behavior of G, one would ide- 
ally want to sample the quasistationary state which is not 
directly available from our simulations for A < Ac- How- 
ever, the finite-time data allow us to extract the short- 
distance behavior up to a crossover distance r^it) given 
by the scaling combination \ii{t/to)r~'^ > 1. Figure 10 
shows the correlation function for infection rates close 
to the critical point. The data show a crossover at a 
length scale of about 200 compatible with the crossover 
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scale identified at the end of Section HID separating 
the pre-asymptotic behavior from the true asymptotic 
long-distance form. Unfortunately, we cannot explore 
the asymptotic region because (i) fluctuations are too 
strong, and (ii) for larger distances, we violate the condi- 
tion \\i{t/t[))r^'^ > 1 discussed above. Thus, the asymp- 
totic functional form of the spatial correlation function 
remains a task for the future. 
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I. Nonequilibrium transition across the lattice 
percolation threshold 

So far, we have considered the generic transition occur- 
ring for impurity concentrations p < Pc- In the present 
subsection, we briefly discuss the nonequilibrium phase 
transition of the contact process across the percolation 
threshold Pc of the underlying lattice. It occurs at suf- 
ficiently large infection rates and is represented by the 
vertical line at pc in the phase diagram shown in Fig. 1. 

Vojta and Lee [41] developed a theory for this tran- 
sition by combining the well-known critical behavior of 
classical percolation with the properties of a supercriti- 
cal contact process on a finite-size cluster. They found 
that the behavior of the contact process on a diluted 
lattice close to the percolation threshold follows the ac- 
tivated scaling scenario of subsection II C. However, the 
critical exponents differ from those of the generic tran- 
sition discussed above; they can be expressed as com- 
binations of the classical lattice percolation exponents 
(3c and which arc known exactly in two space dimen- 
sions. Specifically, the static exponents P and v± are 
identical to their lattice counterparts (3 = Pc = 5/36 
and v±^ = Vc ~ 4/3. The tunneling exponent -0 is 
given by tp = Dc = 2 — pc/vc = 91/48 where Dc is 
the fractal dimension of the critical percolation cluster 
(of the lattice). This implies an extremely small value 
for the exponent controlling the critical density decay, 
P/{iy^i') = 5/91. 

To test this prediction, we performed simulation runs 
at p = Pc and A = 5.0, starting from a fully infected 
lattice. The resulting time evolution of the density of 
active sites (averaged over three samples of linear size 
L = 8000) is presented in Figure 11. The data show the 
behavior predicted in Ref. [41]: a rapid initial density 
decay towards a quasi-stationary density, followed by a 
slow logarithmic time dependence due to the successive 
decays of the contact process on larger and larger perco- 
lation clusters. The asymptotic behavior is compatible 
with the predicted exponent value 6 = 5/91, however, 
determining an accurate value of d from the simulation 
data is impossible because of the extremely slow decay 
and the limited time range we can cover. 



■1 2 3 ,^(,,5 10 

FIG. 11: (Color online) In p(t) versus Inln(f) for p = pc, A = 
5.0, and L — 8000. We have simulated 3 samples, their differ- 
ences are smaller than the line width. The dashed line repre- 
sents the expected logarithmic long-time decay, p ~ [ln(t)]~'' 
with arbitrary prefactor. 



IV. DISCUSSION AND CONCLUSIONS 

In this final section of the paper, we summarize our 
results and discuss some aspects of our data analysis. 
We then briefly consider the diluted contact process in 
three space dimensions. Finally, we relate our findings 
to a recent classification of phase transitions with weak 
disorder. 



A. Summary 

We have performed large-scale Monte-Carlo simula- 
tions of the contact process on a two-dimensional site- 
diluted lattice. We have determined the infection- 
dilution phase diagram und studied both the generic 
nonequilibrium phase transition for dilutions below the 
percolation threshold of the lattice and the transition 
across the lattice percolation threshold. 

Our main results is that the generic transition is con- 
trolled by an infinite-randomness fixed point for all dis- 
order strength investigated, giving rise to ultraslow ac- 
tivated rather than power-law dynamical scaling. Based 
on two types of simulations, (i) spreading of the infec- 
tion from a single seed and (ii) simulations starting from 
a fully infected lattice, we have determined the complete 
critical behavior. Our critical exponent values are sum- 
marized in Table I. Note that the spatial correlation 
length exponent fulfills the Harris criterion [18] inequal- 
ity di'± > 2, as expected in a disordered system [50]. 

We have been able to obtain a rather accurate estimate 
for the finite-size scaling exponent P/v±. However, the 
other exponents are somewhat less accurate despite our 
extensive numerical effort (system sizes up to 8000 x 8000 
sites and times up to 10^°). This is caused by three im- 
portant characteristics of our infinite-randomness critical 
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critical exponent 


value 


value 




(generic) 


(percolation) 




0.96(2) 


5/48 


P 


1.15(15) 


5/36 




1.20(15) 


4/3 




0.51(6) 


91/48 


5 


1.9(2) 


5/91 


e 


0.15(3) 





TABLE I: Critical exponents for the nonequilibrium phase 
transitions of the disordered two-dimensional contact process. 
The values for the generic transition are from this work, the 
values for the percolation transition are from Ref. [41]. The 
numbers in brackets give an error estimate of the last given 
digits. 



point: 

(i) The logarithmic time-dependencies associated with 
activated scahng are not scale-free, they contain a mi- 
croscopic time scale which acts as a strong correction 
to scaling and must be estimated independently (other- 
wise the exponent values will be seriously compromised) . 
Since our simulations cover only a limited range of ln(i), 
our values for to a-rc not very precise which limits the 
accuracy of the exponents 5, Q, and ip. 

(ii) The second problem consists in the extreme nar- 
rowness in A of the asymptotic critical region for the weak 
to moderate disorder strengths we have considered. We 
estimated a relative width of |A — Ad /Ac ~ 10^"^ for di- 
lution p = 0.2. This hampers the analysis of off-critical 
data and thus limits the accuracy of the exponents [3 and 

(iii) The crossover of the time evolution to the asymp- 
totic critical behavior occurs very late. For p = 0.2 the 
crossover time is about 10^. Moreover, the diameter of 
the active cloud spreading from a single seed has reached 
about 200 at that time. Thus, successful simulations re- 
quire both large systems and long times. 

Let us finally comment on the universality of the crit- 
ical exponents. We have studied three different dilution 
values, p = 0.1, 0.2, and 0.3, ranging from weak to moder- 
ate disorder as judged by the shift in the critical infection 
rate with respect to its clean value. While the crossover 
time tx and the microscopic time scale to vary signifi- 
cantly with dilution, the asymptotic critical exponents 
remain constant within their error bars (including the 
rather accurate finite-size scaling exponent (3/i'±). Our 
data thus provide no indication of nonuniversal contin- 
uously varying exponents. However, it must be empha- 
sized that due to the limited precision, some variations 
cannot be rigorously excluded. 

In the case of the transition of the contact process 
across the lattice percolation threshold, our simulation 
data support the theory developed in Ref. [41]: The dy- 
namical critical behavior is activated, with the critical 
exponents given by combinations of the exactly known 
lattice percolation exponents. 




Int 



FIG. 12; (Color online) \nPs{t) and In Ns{t) versus In Int 
for the three-dimensional contact process on a diluted lattice 
{p = 0.3, L = 400, at least 5000 disorder realizations with 128 
trials each). 

B. Three dimensions 

While the existence of infinite-randomness critical 
points in one and two spatial dimensions is well estab- 
lished in the literature, the situation in three dimensions 
is less clear. Within a numerical implementation of the 
strong-disorder renormalization group, Montrunich et al. 
[38] found that the flow is towards strong disorder. How- 
ever, the data were not sufficient for analyzing the critical 
behavior quantitatively. 

We are currently performing Monte-Carlo simulations 
of the three-dimensional contact process on a diluted lat- 
tice analogous to those reported in Sec. III. Preliminary 
results shown in Fig. 12 suggest that the phase tran- 
sition scenario is very similar to that in one and two- 
dimensions: the transition appears to be controlled by 
an infinite-randomness fixed point with activated dynam- 
ical scaling. A detailed quantitative analysis of the criti- 
cal behavior requires significantly longer times and larger 
systems. It will be published elsewhere. 

C. Conclusions 

Recently, a general classification of phase transitions 
in quenched disordered systems with short-range inter- 
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actions has been suggested [39, 51]. It is based on the ef- 
fective dimensionality dcs of the defects (or, equivalently, 
the rare regions.) Three cases can be distinguished. 

(A) If dcff is below the lower critical dimension 
of the problem, the rare region effects are exponentially 
weak, and the critical point is of conventional type. (B) 
In the second class, with doff = d~ , the Griffiths effects 
are of strong power-law type and the critical behavior 
is controlled by an infinite-randomness fixed point with 
activated scaling. (C) For des > d~ , the rare regions 
can undergo the phase transition independently from the 
bulk system. This leads to a destruction of the sharp 
phase transition by smearing [52] . 

The results of this paper agree with this general classi- 
fication scheme as doff = d~ = (this corresponds to rare 
regions being marginal with their life time depending ex- 
ponentially on their size) leading to class B. In contrast, 
the contact process with extended (line or plane) defects 
falls into class C [45, 53]. 

We conclude by pointing out that the unconventional 
behavior found in this paper may explain the striking ab- 
sence of directed percolation scaling [15] in at least some 



of the experiments. However, the extremely slow dynam- 
ics and narrow critical region will prove to be a challenge 
for the verification of the activated scaling scenario not 
just in simulations but also in experiments. We also em- 
phasize that our results are of importance beyond absorb- 
ing state transitions. Since the strong-disorder renormal- 
ization group predicts our transition to be in the univer- 
sality class of the two-dimensional random transverse- 
field Ising model, the critical behavior found here should 
be valid for this entire universality class. 
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